SuperStrict

Framework BRL.Basic

Import BRL.Timer
Import BRL.Random
Import BRL.PNGLoader
Import BRL.Pixmap
Import BRL.StandardIO
Import BRL.GLGraphics
Import Pub.OpenGL
Import Pub.Glew

Const GWIDTH:Int  = 800
Const GHEIGHT:Int = 600

InitGL()

Local SPH1:TSPH     = TSPH.Create( GWIDTH, GHEIGHT, GWIDTH*3.5, GHEIGHT*3.5, 1.0/3.5 )
Local SPH2:TSPH     = TSPH.Create( GWIDTH, GHEIGHT, GWIDTH*5.0, GHEIGHT*5.0, 1.0/5.0 )
Local SPH3:TSPH     = TSPH.Create( GWIDTH, GHEIGHT, GWIDTH*3.0, GHEIGHT*3.0, 1.0/3.0 )

SPH1.SetParticleSize( 9.5 )
SPH2.SetParticleSize( 9.5 )
SPH3.SetParticleSize( 9.5 )

SPH1.SetGravity( 0.0, 0.7 )
SPH2.SetGravity( 0.0, 0.2 )
SPH3.SetGravity( 0.0, 0.4 )

Local Water:TFluid    = TFluid.Create( SPH1 )
Local Solid:TSolid    = TSolid.Create( SPH1 )

Local Mercury:TFluid  = TFluid.Create( SPH2 )
Local GooWater:TFluid = TFluid.Create( SPH2 )
Local Oil:TFluid      = TFluid.Create( SPH2 )

Local Rigid:TRigid    = TRigid.CreateRigidRect( 100, 50, 500, 50, SPH3 )

Rigid.SetColor( 0, 255, 0 )
Solid.SetColor( 255, 0, 0 )
Solid.SetStiffness( 0.9 )

 Mercury.SetColor( 255, 0, 0 )
GooWater.SetColor( 0, 0, 255 )
   Water.SetColor( 0, 0, 255 )
     Oil.SetColor( 0, 255, 0 )

 Mercury.SetRestDensity( 0.5  )
GooWater.SetRestDensity( 0.35 )
   Water.SetRestDensity( 0.35 )
     Oil.SetRestDensity( 0.2  )

 Mercury.SetSurfaceTension( 2.0 )
GooWater.SetSurfaceTension( 2.0 )
   Water.SetSurfaceTension( 0.1 )
     Oil.SetSurfaceTension( 2.0 )

'----------------------- Scenario 1
Local Emitter:TFluidEmitter = TFluidEmitter.Create( Water, GWIDTH*0.1, 230, GWIDTH/6, 5, 3, 0.0, 0.1*1.5 )
Emitter.Trigger()

'Solid.CreateRect( 280, 490, 270, 60 )

CreateBoundaryRect( SPH1, 280, 0, 40, 487 )
CreateBoundaryRect( SPH1, 320, 420, 30, 47 )
CreateBoundaryRect( SPH1, 0, 550, 560, 50 )
CreateBoundaryRect( SPH1, 560, 500, 40, 100 )

'----------------------- Scenario 2
Local Emitter1:TFluidEmitter = TFluidEmitter.Create(  Mercury, GWIDTH*2/5, 230, GWIDTH/6, 5, 3, 0.0, 0.1*3.0 )
Local Emitter2:TFluidEmitter = TFluidEmitter.Create( GooWater, GWIDTH*1/5, 230, GWIDTH/6, 5, 3, 0.0, 0.1*3.0 )
Local Emitter3:TFluidEmitter = TFluidEmitter.Create(      Oil, GWIDTH*3/5, 230, GWIDTH/6, 5, 3, 0.0, 0.1*3.0 )

 Mercury.CreateRect( GWIDTH*0/3, GHEIGHT - 350, GWIDTH/3, 350 )
GooWater.CreateRect( GWIDTH*1/3, GHEIGHT - 350, GWIDTH/3, 350 )
     Oil.CreateRect( GWIDTH*2/3, GHEIGHT - 350, GWIDTH/3, 350 )



Local Scenario:Int = 1

Local Timer:TTimer = CreateTimer( 200 )

While Not ( KeyHit( KEY_ESCAPE ) Or AppTerminate() )
	glClear( GL_COLOR_BUFFER_BIT )
	
	Select Scenario
		Case 1
			If MouseHit( 1 ) Then Solid.CreateRect  (       MouseX(), MouseY(), 50, 150 )
			If MouseHit( 2 ) Then CreateBoundaryRect( SPH1, MouseX(), MouseY(), 50, 50  )
			
			If KeyHit( KEY_SPACE ) Then Emitter.Trigger()
			
			Emitter.Update()
			
			SPH1.Update()
			SPH1.Render()
			
			GLDrawText "Particles: " + SPH1.GetParticleCount(), 100, 0
		Case 2
			If KeyHit( KEY_SPACE ) Then
				Emitter1.Trigger()
				Emitter2.Trigger()
				Emitter3.Trigger()
			EndIf
			
			Emitter1.Update()
			Emitter2.Update()
			Emitter3.Update()
			
			SPH2.Update()
			SPH2.Render()
			
			GLDrawText "Particles: " + SPH2.GetParticleCount(), 100, 0
		Case 3
			If MouseHit( 1 ) Then TRigid.CreateRigidRect( MouseX(), MouseY(), 40, 150, SPH3 )
			
			SPH3.Update()
			SPH3.Render()
			
			GLDrawText "Particles: " + SPH3.GetParticleCount(), 100, 0
	End Select
	
	If KeyHit( KEY_1 ) Then Scenario = 1
	If KeyHit( KEY_2 ) Then Scenario = 2
	If KeyHit( KEY_3 ) Then Scenario = 3
	
	GLDrawText "FPS: " + GetFPS(), 0, 0
	
	Flip 0
	WaitTimer Timer
Wend
End

Global FPSTimer:Int = MilliSecs(), FPSCounter:Int, FPS:Int

Function GetFPS:Int()
	FPSCounter :+ 1
	
	If MilliSecs() - FPSTimer > 1000 Then
		FPS = FPSCounter
		FPSCounter = 0
		
		FPSTimer = MilliSecs()
	EndIf
	
	Return FPS
End Function

Function UserInput( Fluid:TFluid, Scale:Float, SPH:TSPH )
	Local MX:Int = MouseX(), MY:Int = MouseY()
	
	If MouseHit( 1 ) Then Fluid.CreateRect  (      MX*Scale, MY*Scale, 50*Scale, 50*Scale )
	If MouseHit( 2 ) Then CreateBoundaryRect( SPH, MX*Scale, MY*Scale, 50*Scale, 50*Scale )
End Function

Function InitGL()
	GLGraphics GWIDTH, GHEIGHT
	
	glewInit()
	
	glMatrixMode( GL_PROJECTION )
	glLoadIdentity()
	
	glOrtho( 0, GWIDTH, GHEIGHT, 0, -1, 1 )
	
	glMatrixMode( GL_MODELVIEW )
	glLoadIdentity()
	
	glEnable( GL_ALPHA_TEST )
	glAlphaFunc( GL_GREATER, 0.0 )
	
	glClearColor( 0.0, 0.0, 0.0, 1.0 )
End Function

Function CreateBoundaryRect( SPH:TSPH, StartX:Float, StartY:Float, Width:Float, Height:Float )
	StartX :* SPH.InvWorldScale
	StartY :* SPH.InvWorldScale
	Width  :* SPH.InvWorldScale
	Height :* SPH.InvWorldScale
	
	Local X:Float = StartX
	
	While X <= StartX + Width
		Local Y:Float = StartY
		
		While Y <= StartY + Height
			SPH.AddBoundaryParticle( TParticle.CreateBoundaryParticle( SPH, X*TSPH.UNIT_SCALE, Y*TSPH.UNIT_SCALE ) )
			
			Y :+ TSPH.BOUNDARY_SPACING
		Wend
		
		X :+ TSPH.BOUNDARY_SPACING
	Wend
End Function

Type TFluidEmitter
	Field Fluid:TFluid
	
	Field X:Float
	Field Y:Float
	Field Width:Float
	Field Height:Float
	
	Field XSpeed:Float
	Field YSpeed:Float
	
	Field Counter:Int
	Field TimerDelay:Int
	
	Field Active:Int
	
	Method Update()
		If Active Then
			If Counter > TimerDelay Then
				Fluid.CreateRect( X, Y, Width, Height, XSpeed, YSpeed )
				
				Counter = 0
			EndIf
			
			Counter :+ 1
		EndIf
	End Method
	
	Method Trigger()
		Active = Not Active
	End Method
	
	Function Create:TFluidEmitter( Fluid:TFluid, X:Float, Y:Float, W:Float, H:Float, TimerDelay:Int, XSpeed:Float = 0.0, YSpeed:Float = 0.0 )
		Local Emitter:TFluidEmitter = New TFluidEmitter
		
		Emitter.Fluid = Fluid
		Emitter.X = X
		Emitter.Y = Y
		Emitter.Width  = W
		Emitter.Height = H
		Emitter.XSpeed = XSpeed
		Emitter.YSpeed = YSpeed
		Emitter.TimerDelay = TimerDelay
		
		Return Emitter
	End Function
End Type

Type TSPH
	Const     UNIT_SCALE:Float = 0.07
	Const INV_UNIT_SCALE:Float = 1.0/UNIT_SCALE
	
	Const SMOOTHING_LENGTH:Float		= 55*UNIT_SCALE
	Const SMOOTHING_LENGTH_SQ:Float	= SMOOTHING_LENGTH*SMOOTHING_LENGTH
	Const INV_SMOOTHING_LENGTH:Float	= 1.0/SMOOTHING_LENGTH
	
	Const PARTICLE_AREA:Float = ( SMOOTHING_LENGTH*0.3 )^2*Pi
	
	Const    FLUID_SPACING:Float = 85
	Const BOUNDARY_SPACING:Float = 20
	Const    SOLID_SPACING:Float = 25
	Const    RIGID_SPACING:Float = 25
	
	
	Field GWidth:Int
	Field GHeight:Int
	
	Field    WorldScale:Float
	Field InvWorldScale:Float
	
	Field    Timestep:Float
	Field    TimestepSq:Float
	Field InvTimestep:Float
	
	Field ContainerWidth:Float
	Field ContainerHeight:Float
	
	Field    RepulsiveDistance:Float
	Field InvRepulsiveDistance:Float
	Field    RepulsiveDistSq:Float
	Field    RepulsiveD:Float
	
	Field GravityX:Float
	Field GravityY:Float
	
	Field Fluids:TFluid[]
	Field Solids:TSolid[]
	Field Rigids:TRigid[]
	
	Field GridWidth:Int
	Field GridHeight:Int
	
	Field FluidGrid:TParticle[,]
	Field SolidGrid:TParticle[,]
	Field RigidGrid:TParticle[,]
	Field BoundaryGrid:TParticle[,]
	
	Field FluidCount:Int
	Field SolidCount:Int
	Field RigidCount:Int
	Field BoundaryCount:Int
	
	Field FluidParticles:TParticle[]
	Field SolidParticles:TParticle[]
	Field RigidParticles:TParticle[]
	Field BoundaryParticles:TParticle[]
	
	Field ConstraintCount:Int
	Field Constraints:TParticle[][]
	Field ConstraintLength:Float[]
	
	Field ParticleTex:Int
	Field ParticleSize:Float
	
	field tex:TPIxmap
	
	Method New()
		FluidParticles    = New TParticle[ 1000 ]
		SolidParticles    = New TParticle[ 1000 ]
		RigidParticles    = New TParticle[ 1000 ]
		BoundaryParticles = New TParticle[ 1000 ]
		
		Constraints       = [ New TParticle[ 0 ], New TParticle[ 0 ] ]
		Constraints[ 0 ]  = New TParticle[ 1000 ]
		Constraints[ 1 ]  = New TParticle[ 1000 ]
		ConstraintLength  = New Float    [ 1000 ]
		
		SetTimestep( 0.4 )
		SetRepulsion( 30.0, 0.5 )
		SetGravity( 0.0, 0.7 )
		SetParticleSize( 8.0 )
		
		tex = createpixmap(8,8,PF_RGBA8888)
		drawCircleGradient(tex)
		ParticleTex = GLTexFromPixmap(tex) '( LoadPixmapPNG( "Metaball.png" ) )
	End Method
	
	Method SetGraphicsSize( Width:Int, Height:Int )
		GWidth  = Width
		GHeight = Height
	End Method
	
	Method SetContainerSize( Width:Int, Height:Int )
		ContainerWidth  = Width
		ContainerHeight = Height
		
		GridWidth       = Ceil( ContainerWidth *INV_SMOOTHING_LENGTH )
		GridHeight      = Ceil( ContainerHeight*INV_SMOOTHING_LENGTH )
		
		FluidGrid       = New TParticle[ GridWidth, GridHeight ]
		SolidGrid       = New TParticle[ GridWidth, GridHeight ]
		RigidGrid       = New TParticle[ GridWidth, GridHeight ]
		BoundaryGrid    = New TParticle[ GridWidth, GridHeight ]
	End Method
	
	Method SetWorldScale( Scale:Float )
		   WorldScale = Scale
		InvWorldScale = 1.0/Scale
	End Method
	
	Method SetTimestep( NewStep:Float )
		   Timestep   = NewStep
		   TimestepSq = Timestep*Timestep
		InvTimestep   = 1.0/Timestep
	End Method
	
	Method SetRepulsion( Distance:Float, D:Float )
		   RepulsiveDistance = Distance*UNIT_SCALE
		InvRepulsiveDistance = 1.0/RepulsiveDistance
		   RepulsiveDistSq   = RepulsiveDistance*RepulsiveDistance
		
		RepulsiveD:Float = D
	End Method
	
	Method SetGravity( X:Float, Y:Float )
		GravityX = X*UNIT_SCALE
		GravityY = Y*UNIT_SCALE
	End Method
	
	Method SetParticleSize( Size:Float )
		ParticleSize = Size
	End Method
	
	Method Update()
		RedistributeGrid()
		
		Fluid_Reset()
		Solid_Reset()
		Rigid_Reset()
		
		Fluid_UpdateDensity()
		Fluid_CalcForces()
		Solid_CalcForces()
		Rigid_CalcForces()
		
		Constraints_Update()
		
		Fluid_Integrate()
		Solid_Integrate()
		Rigid_Integrate()
	End Method
	
	Method RedistributeGrid()
		For Local X:Int = 0 Until GridWidth
			For Local Y:Int = 0 Until GridHeight
				FluidGrid   [ X, Y ] = Null
				SolidGrid   [ X, Y ] = Null
				RigidGrid   [ X, Y ] = Null
				BoundaryGrid[ X, Y ] = Null
			Next
		Next
		
		MapParticlesToGrid( FluidParticles,    FluidCount,    FluidGrid,    TSPH.INV_SMOOTHING_LENGTH )
		MapParticlesToGrid( SolidParticles,    SolidCount,    SolidGrid,    TSPH.INV_SMOOTHING_LENGTH )
		MapParticlesToGrid( RigidParticles,    RigidCount,    RigidGrid,    TSPH.INV_SMOOTHING_LENGTH )
		MapParticlesToGrid( BoundaryParticles, BoundaryCount, BoundaryGrid, TSPH.INV_SMOOTHING_LENGTH )
	End Method
	
	Method MapParticlesToGrid( Particles:TParticle[], ParticleCount:Int, DestArray:TParticle[,], ScaleFactor:Float )
		For Local I:Int = 0 Until ParticleCount
			Local P:TParticle = Particles[ I ]
			
			Local GridX:Int = Floor( P.PositionX*ScaleFactor )
			Local GridY:Int = Floor( P.PositionY*ScaleFactor )
			
			If GridX >= 0 And GridY >= 0 And GridX < GridWidth And GridY < GridHeight Then
				Local Successor:TParticle = DestArray[ GridX, GridY ]
				
				DestArray[ GridX, GridY ] = P
				
				P.Succ = Successor
			EndIf
		Next
	End Method
	
	Method Fluid_Reset()
		For Local I:Int = 0 Until FluidCount
			Local P:TParticle = FluidParticles[ I ]
			
			P.ForceX = 0.0
			P.ForceY = 0.0
			P.DeltaVelocityX = 0.0
			P.DeltaVelocityY = 0.0
			P.DeltaDensity   = 0.0
		Next
	End Method
	
	Method Fluid_UpdateDensity()
		Local DirX:Int[] = [ -1, -1,  1, 1, -1, 1,  0, 0, 0 ]
		Local DirY:Int[] = [ -1,  1, -1, 1,  0, 0, -1, 1, 0 ]
		
		For Local I:Int = 0 Until FluidCount
			Local P1:TParticle = FluidParticles[ I ]
			
			Local IntX:Int = Int( P1.PositionX*INV_SMOOTHING_LENGTH )
			Local IntY:Int = Int( P1.PositionY*INV_SMOOTHING_LENGTH )
			
			For Local I:Int = 0 To 8
				Local GridX:Int = IntX + DirX[ I ]
				Local GridY:Int = IntY + DirY[ I ]
				
				If GridX < 0 Or GridY < 0 Or GridX >= GridWidth Or GridY >= GridHeight Then Continue
				
				Local P2:TParticle = FluidGrid[ GridX, GridY ]
				
				While P2 <> Null
					If P1.ID > P2.ID 'And P1.Fluid = P2.Fluid Then
						Local DX:Float = ( P1.PositionX - P2.PositionX )
						Local DY:Float = ( P1.PositionY - P2.PositionY )
						
						Local DSQ:Float = DX*DX + DY*DY
						
						If DSQ < SMOOTHING_LENGTH_SQ And DSQ > 0.0 Then
							Local D:Float = Sqr( DSQ )
							
							Local R:Float = SMOOTHING_LENGTH_SQ - DSQ
							
							Local DensityFactor:Float = ( ( P1.VelocityX - P2.VelocityX )*DX + ( P1.VelocityY - P2.VelocityY )*DY )*R*R/D
							
							P1.DeltaDensity :+ P2.Mass*DensityFactor
							P2.DeltaDensity :+ P1.Mass*DensityFactor
						EndIf
					EndIf
					
					P2 = P2.Succ
				Wend
			Next
		Next
		
		For Local I:Int = 0 Until FluidCount
			Local P:TParticle  = FluidParticles[ I ]
			Local Fluid:TFluid = P.Fluid
			
			P.Density :+ P.DeltaDensity*TFluid.D_KERNEL_FACTOR_GRADIENT
			
			If P.Density < P.Fluid.MinDensity Or P.Density > P.Fluid.MaxDensity Then P.Density = Fluid.RestDensity*1.05'; Print P.ID 'CANCER. ELIMINATE.
			
			P.InvDensity   = 1.0/P.Density
			P.InvDensitySq = P.InvDensity*P.InvDensity
			
			P.Mass    = P.Density*PARTICLE_AREA
			P.InvMass = 1.0/P.Mass
			
			P.Pressure = Fluid.P_B*P.Density*( ( P.Density*Fluid.InvRestDensity )^7.0 - 1.0 )
			
			P.SoundSpeed = Sqr( Abs( TFluid.HEAT_RATIO*P.Pressure*P.InvDensity ) )
		Next
	End Method
	
	Method Fluid_CalcForces()
		Local DirX:Int[] = [ -1, -1,  1, 1, -1, 1,  0, 0, 0 ]
		Local DirY:Int[] = [ -1,  1, -1, 1,  0, 0, -1, 1, 0 ]
		
		For Local I:Int = 0 Until FluidCount
			Local P1:TParticle = FluidParticles[ I ]
			
			Local IntX:Int = Int( P1.PositionX*INV_SMOOTHING_LENGTH )
			Local IntY:Int = Int( P1.PositionY*INV_SMOOTHING_LENGTH )
			
			
			For Local I:Int = 0 To 8
				Local GridX:Int = IntX + DirX[ I ]
				Local GridY:Int = IntY + DirY[ I ]
				
				If GridX < 0 Or GridY < 0 Or GridX >= GridWidth Or GridY >= GridHeight Then Continue
				
				Local P2:TParticle = FluidGrid[ GridX, GridY ]
				
				While P2 <> Null
					If P1.ID > P2.ID Then
						Local DX:Float = ( P1.PositionX - P2.PositionX )
						Local DY:Float = ( P1.PositionY - P2.PositionY )
						
						Local DSQ:Float = DX*DX + DY*DY
						
						If DSQ < SMOOTHING_LENGTH_SQ And DSQ > 0.0 Then
							Local D:Float   = Sqr( DSQ )
							
							Local R:Float   = SMOOTHING_LENGTH    - D
							Local R2:Float  = SMOOTHING_LENGTH_SQ - DSQ
							
							Local VDX:Float = P2.VelocityX - P1.VelocityX
							Local VDY:Float = P2.VelocityY - P1.VelocityY
							
							Local PolyKernel:Float = TFluid.D_KERNEL_FACTOR*R2*R2*R2
							
							If P1.Fluid = P2.Fluid Then
								Local TensionForce:Float = -P1.Fluid.SurfaceTension*PolyKernel
								
								Local Tension1:Float = TensionForce*P2.Mass
								Local Tension2:Float = TensionForce*P1.Mass
								
								P1.ForceX :+ DX*Tension1
								P1.ForceY :+ DY*Tension1
								P2.ForceX :- DX*Tension2
								P2.ForceY :- DY*Tension2
							EndIf
							
							Local PressureForce:Float = P1.Pressure*P1.InvDensitySq + P2.Pressure*P2.InvDensitySq
							
							Local InvDensityAverage:Float = 2.0/( P1.Density + P2.Density )
							
							Local DotP:Float = -DX*VDX - DY*VDY
							
							If DotP < 0.0 Then
								Local M:Float = SMOOTHING_LENGTH*DotP/( DSQ + 20.0 )
								
								PressureForce :+ ( -P2.Fluid.V_Alpha*M*( P1.SoundSpeed + P2.SoundSpeed )*0.5 + P2.Fluid.V_Beta*M*M )*InvDensityAverage
							EndIf
														
							PressureForce :* P2.Mass*TFluid.P_KERNEL_FACTOR*R*R/D
							
							DX :* PressureForce
							DY :* PressureForce
							
							P1.ForceX :- DX'*P2.Mass
							P1.ForceY :- DY'*P2.Mass
							P2.ForceX :+ DX'*P1.Mass
							P2.ForceY :+ DY'*P1.Mass
							
							Local KernelFactor:Float = TFluid.XSPH_EPSILON*PolyKernel*InvDensityAverage
							
							VDX :* KernelFactor
							VDY :* KernelFactor
							
							P1.DeltaVelocityX :+ P2.Mass*VDX
							P1.DeltaVelocityY :+ P2.Mass*VDY
							P2.DeltaVelocityX :- P1.Mass*VDX
							P2.DeltaVelocityY :- P1.Mass*VDY
						EndIf
					EndIf
					
					P2 = P2.Succ
				Wend
				
				P2 = BoundaryGrid[ GridX, GridY ]
				
				While P2 <> Null
					Local DX:Float = ( P1.PositionX - P2.PositionX )
					Local DY:Float = ( P1.PositionY - P2.PositionY )
					
					Local DSQ:Float = DX*DX + DY*DY
					
					If DSQ < RepulsiveDistSq Then
						Local InvDSQ:Float   = 1.0/DSQ
						
						Local RSQ:Float = RepulsiveDistSq*InvDSQ
						
						Local RepulsiveFactor:Float = RepulsiveD*( RSQ*RSQ - RSQ )*InvDSQ
						
						P1.ForceX :+ DX*RepulsiveFactor
						P1.ForceY :+ DY*RepulsiveFactor
					EndIf
					
					P2 = P2.Succ
				Wend
				
				P2 = SolidGrid[ GridX, GridY ]
				
				While P2 <> Null
					Local DX:Float = ( P1.PositionX - P2.PositionX )
					Local DY:Float = ( P1.PositionY - P2.PositionY )
					
					Local DSQ:Float = DX*DX + DY*DY
					
					If DSQ < RepulsiveDistSq Then
						Local InvDSQ:Float   = 1.0/DSQ
						
						Local RSQ:Float = RepulsiveDistSq*InvDSQ
						
						Local RepulsiveFactor:Float = RepulsiveD*( RSQ*RSQ - RSQ )*InvDSQ
						
						P1.ForceX :+ DX*RepulsiveFactor*0.5
						P1.ForceY :+ DY*RepulsiveFactor*0.5
						P2.ForceX :- DX*RepulsiveFactor*0.5
						P2.ForceY :- DY*RepulsiveFactor*0.5
					EndIf
					
					P2 = P2.Succ
				Wend
			Next
			
			Local Dist1:Float = Max( ContainerWidth  - P1.PositionX, 0.0 )
			Local Dist2:Float = Max( ContainerHeight - P1.PositionY, 0.0 )
			Local Dist3:Float = Max(                   P1.PositionX, 0.0 )
			Local Dist4:Float = Max(                   P1.PositionY, 0.0 )
			
			'If Dist1 < RepulsiveDistance Then P1.ForceX :- 1.0*( 1.0 - Sqr( Dist1*InvRepulsiveDistance ) )
			'If Dist2 < RepulsiveDistance Then P1.ForceY :- 1.0*( 1.0 - Sqr( Dist2*InvRepulsiveDistance ) )
			'If Dist3 < RepulsiveDistance Then P1.ForceX :+ 1.0*( 1.0 - Sqr( Dist3*InvRepulsiveDistance ) )
			'If Dist4 < RepulsiveDistance Then P1.ForceY :+ 1.0*( 1.0 - Sqr( Dist4*InvRepulsiveDistance ) )
			
			If Dist1 < RepulsiveDistance Then P1.DeltaVelocityX :- 1.0*( 1.0 - Sqr( Dist1*InvRepulsiveDistance ) )
			If Dist2 < RepulsiveDistance Then P1.DeltaVelocityY :- 1.0*( 1.0 - Sqr( Dist2*InvRepulsiveDistance ) )
			If Dist3 < RepulsiveDistance Then P1.DeltaVelocityX :+ 1.0*( 1.0 - Sqr( Dist3*InvRepulsiveDistance ) )
			If Dist4 < RepulsiveDistance Then P1.DeltaVelocityY :+ 1.0*( 1.0 - Sqr( Dist4*InvRepulsiveDistance ) )
			
			P1.ForceX :+ GravityX*P1.Mass
			P1.ForceY :+ GravityY*P1.Mass
		Next
	End Method
	
	Method Fluid_Integrate()
		For Local I:Int = 0 Until FluidCount
			Local P:TParticle = FluidParticles[ I ]
			
			Local OldX:Float = P.PositionX
			Local OldY:Float = P.PositionY
			
			P.PositionX :+ ( 1.0 - P.Fluid.Damping )*( P.PositionX - P.OldX ) + P.DeltaVelocityX*Timestep + TimestepSq*P.ForceX*P.InvMass
			P.PositionY :+ ( 1.0 - P.Fluid.Damping )*( P.PositionY - P.OldY ) + P.DeltaVelocityY*Timestep + TimestepSq*P.ForceY*P.InvMass
			
			P.PositionX = Max( Min( P.PositionX, ContainerWidth  ), 0.0 )
			P.PositionY = Max( Min( P.PositionY, ContainerHeight ), 0.0 )
			
			P.OldX = OldX
			P.OldY = OldY
			
			P.VelocityX = ( P.PositionX - OldX )*InvTimestep
			P.VelocityY = ( P.PositionY - OldY )*InvTimestep
			
			Local DSQ:Float = P.VelocityX*P.VelocityX + P.VelocityY*P.VelocityY
			
			If DSQ > 4.0 Then
				Local Factor:Float = Sqr( 4.0/DSQ )
				
				P.PositionX = OldX
				P.PositionY = OldY
								
				P.VelocityX :* Factor
				P.VelocityY :* Factor
				
				P.PositionX :+ P.VelocityX*Timestep
				P.PositionY :+ P.VelocityY*Timestep
			EndIf
			
			P.ScreenX = P.PositionX*INV_UNIT_SCALE*WorldScale
			P.ScreenY = P.PositionY*INV_UNIT_SCALE*WorldScale
		Next
	End Method
	
	Method Solid_Reset()
		For Local I:Int = 0 Until SolidCount
			Local P:TParticle = SolidParticles[ I ]
			
			P.ForceX = GravityX
			P.ForceY = GravityY
		Next
	End Method
	
	Method Solid_CalcForces()
		Local DirX:Int[] = [ -1, -1,  1, 1, -1, 1,  0, 0, 0 ]
		Local DirY:Int[] = [ -1,  1, -1, 1,  0, 0, -1, 1, 0 ]
		
		For Local I:Int = 0 Until SolidCount
			Local P1:TParticle = SolidParticles[ I ]
			
			Local IntX:Int = Int( P1.PositionX*INV_SMOOTHING_LENGTH )
			Local IntY:Int = Int( P1.PositionY*INV_SMOOTHING_LENGTH )
			
			
			For Local I:Int = 0 To 8
				Local GridX:Int = IntX + DirX[ I ]
				Local GridY:Int = IntY + DirY[ I ]
				
				If GridX < 0 Or GridY < 0 Or GridX >= GridWidth Or GridY >= GridHeight Then Continue
				
				Local P2:TParticle = BoundaryGrid[ GridX, GridY ]
				
				While P2 <> Null
					Local DX:Float = ( P1.PositionX - P2.PositionX )
					Local DY:Float = ( P1.PositionY - P2.PositionY )
					
					Local DSQ:Float = DX*DX + DY*DY
					
					If DSQ < RepulsiveDistSq Then
						Local InvDSQ:Float   = 1.0/DSQ
						
						Local RSQ:Float = RepulsiveDistSq*InvDSQ
						
						Local RepulsiveFactor:Float = RepulsiveD*( RSQ*RSQ - RSQ )*InvDSQ
						
						P1.ForceX :+ DX*RepulsiveFactor
						P1.ForceY :+ DY*RepulsiveFactor
					EndIf
					
					P2 = P2.Succ
				Wend
				
				P2 = SolidGrid[ GridX, GridY ]
				
				While P2 <> Null
					Local DX:Float = ( P1.PositionX - P2.PositionX )
					Local DY:Float = ( P1.PositionY - P2.PositionY )
					
					Local DSQ:Float = DX*DX + DY*DY
					
					If DSQ < RepulsiveDistSq*0.8 And P1.ID > P2.ID Then
						Local InvDSQ:Float   = 1.0/DSQ
						
						Local RSQ:Float = RepulsiveDistSq*0.8*InvDSQ
						
						Local RepulsiveFactor:Float = RepulsiveD*( RSQ*RSQ - RSQ )*InvDSQ
						
						P1.ForceX :+ DX*RepulsiveFactor*0.5
						P1.ForceY :+ DY*RepulsiveFactor*0.5
						P2.ForceX :- DX*RepulsiveFactor*0.5
						P2.ForceY :- DY*RepulsiveFactor*0.5
					EndIf
					
					P2 = P2.Succ
				Wend
			Next
			
			Local Dist1:Float = Max( ContainerWidth  - P1.PositionX, 0.0 )
			Local Dist2:Float = Max( ContainerHeight - P1.PositionY, 0.0 )
			Local Dist3:Float = Max(                   P1.PositionX, 0.0 )
			Local Dist4:Float = Max(                   P1.PositionY, 0.0 )
			
			If Dist1 < RepulsiveDistance Then P1.ForceX :- 1.0*( 1.0 - Sqr( Dist1*InvRepulsiveDistance ) )
			If Dist2 < RepulsiveDistance Then P1.ForceY :- 1.0*( 1.0 - Sqr( Dist2*InvRepulsiveDistance ) )
			If Dist3 < RepulsiveDistance Then P1.ForceX :+ 1.0*( 1.0 - Sqr( Dist3*InvRepulsiveDistance ) )
			If Dist4 < RepulsiveDistance Then P1.ForceY :+ 1.0*( 1.0 - Sqr( Dist4*InvRepulsiveDistance ) )
		Next
	End Method
	
	Method Solid_Integrate()
		For Local I:Int = 0 Until SolidCount
			Local P:TParticle = SolidParticles[ I ]
			
			Local OldX:Float = P.PositionX
			Local OldY:Float = P.PositionY
			
			P.PositionX :+ P.PositionX - P.OldX + TimestepSq*P.ForceX
			P.PositionY :+ P.PositionY - P.OldY + TimestepSq*P.ForceY
			
			P.PositionX = Max( Min( P.PositionX, ContainerWidth  ), 0.0 )
			P.PositionY = Max( Min( P.PositionY, ContainerHeight ), 0.0 )
			
			P.OldX = OldX
			P.OldY = OldY
			
			Local VX:Float = ( P.PositionX - OldX )*InvTimestep
			Local VY:Float = ( P.PositionY - OldY )*InvTimestep
			
			Local DSQ:Float = VX*VX + VY*VY
			
			If DSQ > 4.0 Then
				Local Factor:Float = Sqr( 2.0/DSQ )
				
				P.PositionX = OldX
				P.PositionY = OldY
				
				P.PositionX :+ VX*Factor*Timestep
				P.PositionY :+ VY*Factor*Timestep
			EndIf
			
			P.ScreenX = P.PositionX*INV_UNIT_SCALE*WorldScale
			P.ScreenY = P.PositionY*INV_UNIT_SCALE*WorldScale
		Next
	End Method
	
	Method Rigid_Reset()
		For Local I:Int = 0 Until RigidCount
			Local P:TParticle = RigidParticles[ I ]
			
			P.ForceX = 0.0
			P.ForceY = 0.0
		Next
	End Method
	
	Method Rigid_CalcForces()
		Local DirX:Int[] = [ -1, -1,  1, 1, -1, 1,  0, 0, 0 ]
		Local DirY:Int[] = [ -1,  1, -1, 1,  0, 0, -1, 1, 0 ]
		
		For Local I:Int = 0 Until RigidCount
			Local P:TParticle = RigidParticles[ I ]
			
			Local Dist1:Float = Max( ContainerWidth  - P.PositionX, 0.0 )
			Local Dist2:Float = Max( ContainerHeight - P.PositionY, 0.0 )
			Local Dist3:Float = Max(                   P.PositionX, 0.0 )
			Local Dist4:Float = Max(                   P.PositionY, 0.0 )
			
			If Dist1 <= RepulsiveDistance Then P.ForceX :- P.Rigid.Repulsion*( 1.0 - Sqr( Dist1*InvRepulsiveDistance ) )
			If Dist2 <= RepulsiveDistance Then P.ForceY :- P.Rigid.Repulsion*( 1.0 - Sqr( Dist2*InvRepulsiveDistance ) )
			If Dist3 <= RepulsiveDistance Then P.ForceX :+ P.Rigid.Repulsion*( 1.0 - Sqr( Dist3*InvRepulsiveDistance ) )
			If Dist4 <= RepulsiveDistance Then P.ForceY :+ P.Rigid.Repulsion*( 1.0 - Sqr( Dist4*InvRepulsiveDistance ) )
			
			If P.PositionX > ContainerWidth  Then Rigid_ResolvePenetration( P, -1,  0,  P.PositionX - ContainerWidth  )
			If P.PositionY > ContainerHeight Then Rigid_ResolvePenetration( P,  0, -1,  P.PositionY - ContainerHeight )
			If P.PositionX < 0               Then Rigid_ResolvePenetration( P,  1,  0, -P.PositionX )
			If P.PositionY < 0               Then Rigid_ResolvePenetration( P,  0,  1, -P.PositionY )
		Next
	End Method
	
	Method Rigid_Integrate()
		For Local Rigid:TRigid = EachIn Rigids
			Rigid.Integrate()
		Next
	End Method
	
	Method Rigid_ResolvePenetration( P:TParticle, NormalX:Float, NormalY:Float, Depth:Float, E:Float = 0.0 )
		Local RX:Float = P.PositionY - P.Rigid.CenterY
		Local RY:Float = P.Rigid.CenterX - P.PositionX
		
		Local VelocityX:Float = P.Rigid.VelocityX + RX*P.Rigid.AngularSpeed
		Local VelocityY:Float = P.Rigid.VelocityY + RY*P.Rigid.AngularSpeed
		
		Local NormalVelocity:Float = VelocityX*NormalX + VelocityY*NormalY
		
		If NormalVelocity < 0.0 Then
			Local PerpProduct:Float = RX*NormalX + RY*NormalY
			
			Local J:Float = -( 1 + E )*NormalVelocity/( ( NormalX*NormalX + NormalY*NormalY )*P.Rigid.InvMass + PerpProduct*PerpProduct*P.Rigid.InvInertia )
			
			P.Rigid.VelocityX    :+ NormalX*J*P.Rigid.InvMass
			P.Rigid.VelocityY    :+ NormalY*J*P.Rigid.InvMass
			
			P.Rigid.AngularSpeed :+ PerpProduct*J*P.Rigid.InvInertia
		EndIf
	End Method
	
	Method Constraints_Update()
		For Local I:Int = 0 Until ConstraintCount
			Local P1:TParticle = Constraints[ 0 ][ I ]
			Local P2:TParticle = Constraints[ 1 ][ I ]
			
			Local dX:Float = P2.PositionX - P1.PositionX
			Local dY:Float = P2.PositionY - P1.PositionY
			
			Local dLength:Float = Sqr( dX*dX + dY*dY )
			Local Diff:Float    = ( dLength - ConstraintLength[ I ] )/dLength
			
			Local Stiffness:Float = P1.Solid.Stiffness
			
			P1.PositionX :+ dX*Diff*0.5*Stiffness
			P1.PositionY :+ dY*Diff*0.5*Stiffness
			P2.PositionX :- dX*Diff*0.5*Stiffness
			P2.PositionY :- dY*Diff*0.5*Stiffness
		Next
	End Method
	
	Method Render()
		glEnable( GL_TEXTURE_2D )
		glBindTexture( GL_TEXTURE_2D, ParticleTex )
		
		glEnable( GL_BLEND )
		glBlendFunc( GL_ONE, GL_ONE )
		
		glBegin( GL_QUADS )
			RenderParticleRects( FluidParticles,    FluidCount,    ParticleSize )
		glEnd()
		
		glBegin( GL_QUADS )
			RenderParticleRects( SolidParticles,    SolidCount,    ParticleSize )
			RenderParticleRects( RigidParticles,    RigidCount,    ParticleSize )
		glEnd()
		
		glDisable( GL_TEXTURE_2D )
		
		glBlendFunc( GL_ZERO, GL_DST_COLOR )
		glColor3f( 1.0, 1.0, 1.0 )
		
		glBegin( GL_QUADS )
			For Local I:Int = 0 To 1
				glVertex2f( 0.0, 0.0 )
				glVertex2f( GWidth, 0.0 )
				glVertex2f( GWidth, GHeight )
				glVertex2f( 0.0, GHeight )
			Next
		glEnd()
		
		'glColor3f( 1.0, 0.5, 0.5 )
		'glBlendFunc( GL_DST_COLOR, GL_ZERO )
		
		'glBegin( GL_QUADS )
		'	glVertex2f( 0, 0 )
		'	glVertex2f( GWidth, 0 )
		'	glVertex2f( GWidth, GHeight )
		'	glVertex2f( 0, GHeight )
		'glEnd()
		
		glDisable( GL_BLEND )
		
		glBegin( GL_QUADS )
			RenderParticleRects( BoundaryParticles, BoundaryCount, ParticleSize )
		glEnd()
		
		glColor3f( 1.0, 1.0, 1.0 )
	End Method
	
	Method RenderParticleRects( Particles:TParticle[], ParticleCount:Int, Size:Float )
		glColor3f( 1.0, 1.0, 1.0 )
		
		For Local I:Int = 0 Until ParticleCount
			Local P:TParticle = Particles[ I ]
			
			glColor3f( P.ColorR*0.00390625, P.ColorG*0.00390625, P.ColorB*0.00390625 )
			
			glTexCoord2f( 0.0, 0.0 ); glVertex2f( P.ScreenX - Size, P.ScreenY - Size )
			glTexCoord2f( 1.0, 0.0 ); glVertex2f( P.ScreenX + Size, P.ScreenY - Size )
			glTexCoord2f( 1.0, 1.0 ); glVertex2f( P.ScreenX + Size, P.ScreenY + Size )
			glTexCoord2f( 0.0, 1.0 ); glVertex2f( P.ScreenX - Size, P.ScreenY + Size )
		Next
	End Method
	
	Method GetParticleCount:Int()
		Return FluidCount + SolidCount + RigidCount + BoundaryCount
	End Method
	
	Method AddFluid( Fluid:TFluid )
		Fluids :+ [ Fluid ]
	End Method
	
	Method AddSolid( Solid:TSolid )
		Solids :+ [ Solid ]
	End Method
	
	Method AddRigid( Rigid:TRigid )
		Rigids :+ [ Rigid ]
	End Method
	
	Method AddFluidParticle( P:TParticle )
		If FluidCount = FluidParticles.Length Then FluidParticles = FluidParticles[ .. FluidParticles.Length*2 ]
		
		FluidParticles[ FluidCount ] = P
		FluidCount :+ 1
	End Method
	
	Method AddSolidParticle( P:TParticle )
		If SolidCount = SolidParticles.Length Then SolidParticles = SolidParticles[ .. SolidParticles.Length*2 ]
		
		SolidParticles[ SolidCount ] = P
		SolidCount :+ 1
	End Method
	
	Method AddRigidParticle( P:TParticle )
		If RigidCount = RigidParticles.Length Then RigidParticles = RigidParticles[ .. RigidParticles.Length*2 ]
		
		RigidParticles[ RigidCount ] = P
		RigidCount :+ 1
	End Method
	
	Method AddBoundaryParticle( P:TParticle )
		If BoundaryCount = BoundaryParticles.Length Then BoundaryParticles = BoundaryParticles[ .. BoundaryParticles.Length*2 ]
		
		BoundaryParticles[ BoundaryCount ] = P
		BoundaryCount :+ 1
	End Method
	
	Method AddConstraint( P1:TParticle, P2:TParticle )
		If ConstraintCount = ConstraintLength.Length Then
			Constraints[ 0 ] = Constraints[ 0 ][ .. ConstraintLength.Length*2 ]
			Constraints[ 1 ] = Constraints[ 1 ][ .. ConstraintLength.Length*2 ]
			ConstraintLength = ConstraintLength[ .. ConstraintLength.Length*2 ]
		EndIf
		
		Constraints[ 0 ][ ConstraintCount ] = P1
		Constraints[ 1 ][ ConstraintCount ] = P2
		ConstraintLength[ ConstraintCount ] = Sqr( ( P1.PositionX - P2.PositionX )*( P1.PositionX - P2.PositionX ) + ..
		                                           ( P1.PositionY - P2.PositionY )*( P1.PositionY - P2.PositionY ) )
		
		ConstraintCount :+ 1
	End Method
	
	Function Create:TSPH( GWidth:Int, GHeight:Int, FluidWidth:Int, FluidHeight:Int, FluidScale:Float )
		Local SPH:TSPH = New TSPH
		
		SPH.SetGraphicsSize( GWidth, GHeight )
		SPH.SetContainerSize( FluidWidth*UNIT_SCALE, FluidHeight*UNIT_SCALE )
		SPH.SetWorldScale( FluidScale )
		
		Return SPH
	End Function
End Type

Type TSolid
	Field Integrator:TSPH
	
	Field Stiffness:Float
	
	Field ColorR:Byte
	Field ColorG:Byte
	Field ColorB:Byte
	
	Method New()
		SetStiffness( 0.9 )
	End Method
	
	Method SetIntegrator( SPH:TSPH )
		Integrator = SPH
		
		Integrator.AddSolid( Self )
	End Method
	
	Method SetStiffness( Value:Float )
		Stiffness = Value
	End Method
	
	Method SetColor( R:Byte, G:Byte, B:Byte )
		ColorR = R
		ColorG = G
		ColorB = B
	End Method
	
	Method GetColor( R:Byte Var, G:Byte Var, B:Byte Var )
		R = ColorR
		G = ColorG
		B = ColorB
	End Method
	
	Method CreateRect( StartX:Float, StartY:Float, Width:Float, Height:Float, XSpeed:Float = 0.0, YSpeed:Float = 0.0 )
		StartX :* Integrator.InvWorldScale
		StartY :* Integrator.InvWorldScale
		Width  :* Integrator.InvWorldScale
		Height :* Integrator.InvWorldScale
		
		Local Grid:TParticle[ Ceil( Width/TSPH.SOLID_SPACING ), Ceil( Height/TSPH.SOLID_SPACING ) ]
		Local X:Float = StartX
		
		For Local X:Float = 0 Until Width Step TSPH.SOLID_SPACING
			For Local Y:Float = 0 Until Height Step TSPH.SOLID_SPACING
				Local P:TParticle = TParticle.CreateSolidParticle( ( X + StartX )*TSPH.UNIT_SCALE, ..
				                                                   ( Y + StartY )*TSPH.UNIT_SCALE, Self )
				
				P.OldX :- XSpeed
				P.OldY :- YSpeed
				
				Grid[ X/TSPH.SOLID_SPACING, Y/TSPH.SOLID_SPACING ] = P
				
				Integrator.AddSolidParticle( P )
			Next
		Next
		
		For Local X:Int = 0 Until Grid.Dimensions()[ 0 ] - 1
			For Local Y:Int = 0 Until Grid.Dimensions()[ 1 ] - 1
				Integrator.AddConstraint( Grid[ X + 0, Y + 0 ], Grid[ X + 1, Y + 0 ] )
				Integrator.AddConstraint( Grid[ X + 0, Y + 0 ], Grid[ X + 0, Y + 1 ] )
				Integrator.AddConstraint( Grid[ X + 0, Y + 1 ], Grid[ X + 1, Y + 1 ] )
				Integrator.AddConstraint( Grid[ X + 1, Y + 0 ], Grid[ X + 1, Y + 1 ] )
			Next
		Next
		
		For Local X:Int = 0 Until Grid.Dimensions()[ 0 ] - 2
			For Local Y:Int = 0 Until Grid.Dimensions()[ 1 ] - 2
				Integrator.AddConstraint( Grid[ X + 0, Y + 0 ], Grid[ X + 2, Y + 2 ] )
				Integrator.AddConstraint( Grid[ X + 0, Y + 2 ], Grid[ X + 2, Y + 0 ] )
			Next
		Next
	End Method
	
	Function Create:TSolid( SPH:TSPH )
		Local Solid:TSolid = New TSolid
		
		Solid.SetIntegrator( SPH )
		
		Return Solid
	End Function
End Type

Type TFluid
	Const D_KERNEL_FACTOR:Float			= 315.0/( 64.0*Pi*TSPH.SMOOTHING_LENGTH^9 )	'Poly6
	Const P_KERNEL_FACTOR:Float			= -45.0/(      Pi*TSPH.SMOOTHING_LENGTH^6 )	'Spiky (gradient)
	Const D_KERNEL_FACTOR_GRADIENT:Float	= -6.0*D_KERNEL_FACTOR
	
	Const HEAT_RATIO:Float   = 1.0
	
	Const XSPH_EPSILON:Float = 0.5
	
	
	Field Integrator:TSPH
	
	Field Damping:Float
	
	Field SurfaceTension:Float
	
	Field V_Alpha:Float
	Field V_Beta:Float
	
	Field    RestDensity:Float
	Field InvRestDensity:Float
	
	Field MinDensity:Float
	Field MaxDensity:Float
	
	Field P_B:Float
	
	Field ColorR:Byte
	Field ColorG:Byte
	Field ColorB:Byte
	
	Method New()
		SetDamping( 0.0015 )
		SetViscosity( 5.0, 10.0 )
		SetRestDensity( 0.35 )
		SetPressureFactor( 5.0 )
	End Method
	
	Method SetIntegrator( SPH:TSPH )
		Integrator = SPH
		
		Integrator.AddFluid( Self )
	End Method
	
	Method SetDamping( Value:Float )
		Damping = Value
	End Method
	
	Method SetSurfaceTension( Value:Float )
		SurfaceTension = Value
	End Method
	
	Method SetViscosity( Alpha:Float, Beta:Float )
		V_Alpha = Alpha
		V_Beta  = Beta
	End Method
	
	Method SetRestDensity( Density:Float )
		   RestDensity = Density
		InvRestDensity = 1.0/Density
		
		MinDensity = Density*0.85
		MaxDensity = Density*1.15
	End Method
	
	Method SetPressureFactor( Factor:Float )
		P_B = Factor
	End Method
	
	Method SetColor( R:Byte, G:Byte, B:Byte )
		ColorR = R
		ColorG = G
		ColorB = B
	End Method
	
	Method GetColor( R:Byte Var, G:Byte Var, B:Byte Var )
		R = ColorR
		G = ColorG
		B = ColorB
	End Method
	
	Method CreateRect( StartX:Float, StartY:Float, Width:Float, Height:Float, XSpeed:Float = 0.0, YSpeed:Float = 0.0 )
		Local Spacing:Float = TSPH.FLUID_SPACING*0.35' + Min( TSPH.FLUID_SPACING*( RestDensity - 0.4 )*0.3, 0.0 )
		
		StartX :* Integrator.InvWorldScale
		StartY :* Integrator.InvWorldScale
		Width  :* Integrator.InvWorldScale
		Height :* Integrator.InvWorldScale
		
		Local X:Float = StartX
		
		While X <= StartX + Width
			Local Y:Float = StartY
			
			While Y <= StartY + Height
				Local P:TParticle = TParticle.CreateFluidParticle( X*TSPH.UNIT_SCALE, Y*TSPH.UNIT_SCALE, Self )
				
				P.OldX :- XSpeed
				P.OldY :- YSpeed
				
				Integrator.AddFluidParticle( P )
				
				Y :+ Spacing
			Wend
			
			X :+ Spacing
		Wend
	End Method
	
	Function Create:TFluid( SPH:TSPH )
		Local Fluid:TFluid = New TFluid
		
		Fluid.SetIntegrator( SPH )
		
		Return Fluid
	End Function
End Type

Type TRigid
	Const MIN_VERTICES:Int = 64
	
	Const RAD_TO_DEG:Float = 360.0/( 2.0*Pi )
	Const DEG_TO_RAD:Float = 2.0*Pi/360.0
	
	Global CurrBody:TRigid
	
	Field Integrator:TSPH
	
	Field ColorR:Byte
	Field ColorG:Byte
	Field ColorB:Byte
	
	Field CenterX:Float
	Field CenterY:Float
	
	Field Repulsion:Float
	
	Field    Mass:Float
	Field InvMass:Float
	
	Field    Inertia:Float
	Field InvInertia:Float
	
	Field VelocityX:Float
	Field VelocityY:Float
	
	Field AccelerationX:Float
	Field AccelerationY:Float
	
	Field Angle:Float
	Field AngularSpeed:Float
	
	Field VertexCount:Int
	Field Vertices:TParticle[]
	
	Method New()
		Vertices = New TParticle[ MIN_VERTICES ]
		
		SetColor( 0, 255, 0 )
	End Method
	
	Method SetIntegrator( SPH:TSPH )
		Integrator = SPH
		
		Integrator.AddRigid( Self )
	End Method
	
	Method SetRepulsion( Value:Float )
		Repulsion = Value
	End Method
	
	Method SetColor( R:Byte, G:Byte, B:Byte )
		ColorR = R
		ColorG = G
		ColorB = B
	End Method
	
	Method GetColor( R:Byte Var, G:Byte Var, B:Byte Var )
		R = ColorR
		G = ColorG
		B = ColorB
	End Method
	
	Method AddParticle( P:TParticle )
		If Vertices.Length = VertexCount Then Vertices = Vertices[ .. VertexCount*2 ]
		
		Vertices[ VertexCount ] = P
		
		VertexCount :+ 1
	End Method
	
	Method Init()
		For Local I:Int = 0 Until VertexCount
			Local P:TParticle = Vertices[ I ]
			
			Mass :+ P.Mass
			
			CenterX :+ P.PositionX*P.Mass
			CenterY :+ P.PositionY*P.Mass
		Next
		
		InvMass = 1.0/Mass
		
		CenterX :* InvMass
		CenterY :* InvMass
		
		For Local I:Int = 0 Until VertexCount
			Local P:TParticle = Vertices[ I ]
			
			Local DistX:Float = P.PositionX - CenterX
			Local DistY:Float = P.PositionY - CenterY
			
			Local DistSq:Float = DistX*DistX + DistY*DistY
			
			Inertia :+ P.Mass*DistSq
			
			P.Radius = Sqr( DistSq )
			P.Angle  = ATan2( DistY, DistX )
			
			P.LocalX = P.PositionX - CenterX
			P.LocalY = P.PositionY - CenterY
		Next
		
		InvInertia = 1.0/Inertia
		
		SetRepulsion( Mass*0.0008316*1.0 )
	End Method
	
	Method Integrate()
		Local Torque:Float
		
		Local ForceX:Float = Integrator.GravityX*Mass
		Local ForceY:Float = Integrator.GravityY*Mass
		
		For Local I:Int = 0 Until VertexCount
			Local P:TParticle = Vertices[ I ]
			
			Local DiffX:Float = CenterX - P.PositionX
			Local DiffY:Float = CenterY - P.PositionY
			
			Local Tau:Float = DiffX*P.ForceY - DiffY*P.ForceX
			
			Torque :+ Tau/Inertia
			
			ForceX :+ P.ForceX
			ForceY :+ P.ForceY
		Next
		
		VelocityX    :+ Integrator.Timestep*ForceX/Mass
		VelocityY    :+ Integrator.Timestep*ForceY/Mass
		
		CenterX      :+ VelocityX*Integrator.Timestep
		CenterY      :+ VelocityY*Integrator.Timestep
		
		AngularSpeed :+ Torque*Integrator.Timestep
		Angle        :+ RAD_TO_DEG*AngularSpeed*Integrator.Timestep
		
		For Local I:Int = 0 Until VertexCount
			Local P:TParticle = Vertices[ I ]
			
			P.PositionX =  Cos( Angle + P.Angle )*P.Radius + CenterX
			P.PositionY = -Sin( Angle + P.Angle )*P.Radius + CenterY
			
			P.ScreenX = P.PositionX*TSPH.INV_UNIT_SCALE*Integrator.WorldScale
			P.ScreenY = P.PositionY*TSPH.INV_UNIT_SCALE*Integrator.WorldScale
		Next
	End Method
	
	Function BeginBody( SPH:TSPH )
		CurrBody = New TRigid
		
		CurrBody.SetIntegrator( SPH )
	End Function
	
	Function AddVertex:TParticle( X:Float, Y:Float, Mass:Float )
		Local P:TParticle = TParticle.CreateRigidParticle( X, Y, Mass, CurrBody )
		
		CurrBody.AddParticle( P )
		
		Return P
	End Function
	
	Function EndBody:TRigid()
		CurrBody.Init()
		
		CurrBody.Integrator.AddRigid( CurrBody )
		
		Return CurrBody
	End Function
	
	Function CreateRigidRect:TRigid( StartX:Float, StartY:Float, Width:Float, Height:Float, SPH:TSPH )
		StartX :* SPH.InvWorldScale
		StartY :* SPH.InvWorldScale
		Width  :* SPH.InvWorldScale
		Height :* SPH.InvWorldScale
		
		Local X:Float = StartX
		
		BeginBody( SPH )
		
		While X <= StartX + Width
			Local Y:Float = StartY
			
			While Y <= StartY + Height
				Local P:TParticle = AddVertex( X*TSPH.UNIT_SCALE, Y*TSPH.UNIT_SCALE, 1.0 )
				
				SPH.AddRigidParticle( P )
				
				Y :+ TSPH.RIGID_SPACING
			Wend
			
			X :+ TSPH.RIGID_SPACING
		Wend
		
		Return EndBody()
	End Function
End Type

Type TParticle
	Global Count:Int
	
	Global FluidCount:Int
	Global SolidCount:Int
	Global RigidCount:Int
	Global BoundaryCount:Int
	
	Field Fluid:TFluid
	Field Solid:TSolid
	Field Rigid:TRigid
	
	Field ColorR:Byte
	Field ColorG:Byte
	Field ColorB:Byte
	
	Field ScreenX:Float
	Field ScreenY:Float
	
	Field PositionX:Float
	Field PositionY:Float
	
	Field ForceX:Float
	Field ForceY:Float
	
	Field Mass:Float
	Field InvMass:Float
	
	Field Succ:TParticle
	Field ID:Int
	
	'Fluid
	Field OldX:Float
	Field OldY:Float
	
	Field VelocityX:Float
	Field VelocityY:Float
	
	Field DeltaVelocityX:Float
	Field DeltaVelocityY:Float
	
	Field Density:Float
	Field DeltaDensity:Float
	Field InvDensity:Float
	Field InvDensitySq:Float
	
	Field Pressure:Float
	Field SoundSpeed:Float
	
	'Rigid
	Field Radius:Float
	Field Angle:Float
	
	Field LocalX:Float
	Field LocalY:Float
	
	
	Function CreateFluidParticle:TParticle( X:Float, Y:Float, Fluid:TFluid )
		Local P:TParticle = New TParticle
			P.Fluid = Fluid
			P.PositionX = X
			P.PositionY = Y
			P.OldX = X
			P.OldY = Y
			P.ScreenX = X*TSPH.INV_UNIT_SCALE*Fluid.Integrator.WorldScale
			P.ScreenY = Y*TSPH.INV_UNIT_SCALE*Fluid.Integrator.WorldScale
			P.Density = Fluid.RestDensity
			P.Mass = TSPH.PARTICLE_AREA*Fluid.RestDensity
			P.ID = Count
		
		Fluid.GetColor( P.ColorR, P.ColorG, P.ColorB )
		
		Count :+ 1
		FluidCount :+ 1
		
		Return P
	End Function
	
	Function CreateSolidParticle:TParticle( X:Float, Y:Float, Solid:TSolid )
		Local P:TParticle = New TParticle
			P.Solid = Solid
			P.PositionX = X
			P.PositionY = Y
			P.OldX = X
			P.OldY = Y
			P.ScreenX = X*TSPH.INV_UNIT_SCALE*Solid.Integrator.WorldScale
			P.ScreenY = Y*TSPH.INV_UNIT_SCALE*Solid.Integrator.WorldScale
			P.ID = Count
		
		Solid.GetColor( P.ColorR, P.ColorG, P.ColorB )
		
		Count :+ 1
		SolidCount :+ 1
		
		Return P
	End Function
	
	Function CreateRigidParticle:TParticle( X:Float, Y:Float, Mass:Float, Rigid:TRigid )
		Local P:TParticle = New TParticle
			P.Rigid = Rigid
			P.PositionX = X
			P.PositionY = Y
			P.Mass = Mass
			P.ScreenX = X*TSPH.INV_UNIT_SCALE*Rigid.Integrator.WorldScale
			P.ScreenY = Y*TSPH.INV_UNIT_SCALE*Rigid.Integrator.WorldScale
			P.ID = Count
		
		Rigid.GetColor( P.ColorR, P.ColorG, P.ColorB )
		
		Count :+ 1
		RigidCount :+ 1
		
		Return P
	End Function
	
	Function CreateBoundaryParticle:TParticle( SPH:TSPH, X:Float, Y:Float )
		Local P:TParticle = New TParticle
			P.PositionX = X
			P.PositionY = Y
			P.ScreenX = X*TSPH.INV_UNIT_SCALE*SPH.WorldScale
			P.ScreenY = Y*TSPH.INV_UNIT_SCALE*SPH.WorldScale
			P.ID = Count
		
		P.ColorR = 255
		P.ColorG = 128
		P.ColorB = 128
		
		Count :+ 1
		BoundaryCount :+ 1
		
		Return P
	End Function
End Type

Function drawCircleGradient(px:TPixmap)
	Local siz:Float = px.width
	Local center:Float = (siz-1) / 2.0
   	For Local y:Int=0 Until siz
   	For Local x:Int=0 Until siz
     	 	Local t:Float = Sqr( (center-x)*(center-x) + (center-y)*(center-y) ) / center
		'Local c:Int = 255-(t)*255
		Local c:Int = 250-(t^0.6)*255
		If( c < 0 ) c = 0
		If( c > 255 ) c = 255
		
		Local argb:Int = ( c Shl 24) + c Shl 16 + c Shl 8 + c
		WritePixel px,x,y,argb      	
   	Next
	Next 
End Function
